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Abstract 



We consider a family of multi-phase Stefan problems for a certain 1-d model of cell-to-cell 
adhesion and diffusion, which takes the form of a nonlinear forward-backward parabolic equation. 
In each material phase the cell density stays either high or low, and phases are connected by 
jumps across an 'unstable' interval. We develop an existence theory for such problems which 
allows for the annihilation of phases and the subsequent continuation of solutions. Stability 
results for the long-time behaviour of solutions are also obtained, and, where necessary, the 
analysis is complemented by numerical simulations. 

1 Introduction 

In this paper, we give further consideration to the 1-d continuum model for adhesion/diffusion of 
biological cells developed by Anguige and Schmeiser in [3] , which took the form of the nonlinear 
diffusion equation 



for the scaled cell density p(x,t) £ [0,1], and the adhesion coefficient a 6 [0,1], the boundary 
condition being just ||=0ata; = 0,l. 

These equations were obtained as the formal continuum limit of the fundamental biased-random- 
walk model 




(1) 



with quadratic diffusivity 




(2) 



dfh 
dt 

with transitional probabilities 




(3) 




(1 - ft±i)(l - ap iTl )/h 2 , 



(4) 
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on a lattice of points Xi = ih, by taking Taylor expansions about Xi, and letting h — ► 0. In the 
derivation of this equation, h was interpreted as a (microscopic) measure of cell size. 

We recall from [3J that (HI)-© is globally well posed if a < §. If, on the other hand, a > § then 
(fl])-© is ill-posed iff the initial density profile protrudes into the 'unstable' interval 

/„ = (pW(a)) == f 2a ' V f a ' 3) , 2a+ ^ (4a ' 3) l C [1/3,1], (5) 
\ da da I 

since in that case D(p) is positive iff p ^ / Q , and positivity is preserved by the Maximum Principle. 

For completeness, note that in the borderline case a = |, equation (fT]) is just the porous- medium 
equation with quadratic diffusivity and possible change of sign about p = § . For initial data which 
stays away from p = | (either above or below), (fT]) is uniformly parabolic, and global existence of a 
smooth solution follows as for a < |, while for degenerate initial data one is merely guaranteed a 
(unique) globally existing weak solution [5]. 

The ill-posedness of (JXJ) for a > § is related to the presence of fine (wavelength 0(h)) spatial 
oscillations, as well as plateau formation, in solutions of the discrete system ([3]), and the absence 
of a straightforward existence theory for (fT]) leads one to ask just what model should be taken as a 
reasonable continuum limit of ([3j in the high-adhesion regime. 

One approach, and the one we shall adopt in this paper, is to circumvent the problem of ill- 
posedness by simply declaring that p- values in I a are forbidden, and considering solutions to (fl])-© 
which may jump across I a (possibly multiple times), but which are otherwise smooth. Mathemati- 
cally, one is then dealing with a kind of (multi-phase) Stefan problem for the density p(x, t) and the 
jump locations Sj(t), such that the Si are dynamically determined by local conservation of mass, or, 
in other words, by the Rankine-Hugoniot condition. 

In [3] , the analysis of ([3J was aided by considering higher-order modifications of the leading-order 
equation UJ. One such 0(h 2 )-modification takes the form of the fourth-order PDE 

where the cubic K(p) is a primitive for D(p). This equation is rather similar to the viscous Cahn- 
Hilliard equation [3] , and is a regularisation of |T]) in the sense that it is (at least locally) well-posed 
on S 1 , for each fixed value of the microscopic parameter h, and for all a < 1. Presumably, solutions 
continue to exist globally, as for Cahn-Hilliard, but a proof is currently lacking. 

The steady-state equation for ([6]) can, after a change of variables, be written as a Hamiltonian 
dynamical system, and amongst the solutions there is, for each a, a unique heteroclinic cycle. 
These heteroclinic cycles correspond to (two-level) plateau solutions of ([6]), are close (for small h) 
to square-wave weak solutions of (fl]), and their critical points, denoted by pi(a) and P2(a), such 
that pi < p b < < p2 and K(pi) = K(p2), match very well the numerically-observed long-time 
plateau values in solutions of ([3]) (see [3j). For these reasons, we demand in our Stefan-problem 
framework that any jumps across I a should connect p\(a) to P2{oc)- Furthermore, in order to avoid 
the degeneracies at dl a , we require that the initial data satisfy p < p b in low-density phases, and 
p > p$ in high-density ones. 

We emphasise that p\ and p2 are determined by the particular choice of microscopic model 
©-(H]); other model choices are possible, and these will result in different p- values. 

The paper is organised as follows. In Section 2, we develop a partial existence theory and 
perform a steady-state analysis for the simplest Stefan problem, namely, that for which there is only 
a single discontinuity in the density; solution behaviour is further clarified with the aid of several 
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Figure 1: A typical density profile for the 1-jump Stefan problem, SP\. 



numerical simulations. In Section 3, we extend the analysis to the general multi-phase case, which, 
in particular, allows for the annihilation of phases via coalescence events. Finally, in the Appendices, 
we collect a number of results from classical parabolic theory which are used throughout the paper. 

2 The 1-jump problem 

We begin the analysis by considering the simplest possible case, in which there are just two phases, 
connected by a single jump from pi to pn\ this situation is depicted in Figure 1. For defmiteness, 
we will assume that the low-density phase lies to the left, and the high-density phase to the right, 
of the discontinuity; the converse arrangement can, of course, also be treated. 
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2.1 Formulation 



Let the location of the jump discontinuity be denoted by s(t). The Stefan problem for a given a > 
then consists of looking for a function p(x, t) on [0, 1] x [0, T] which satisfies 



subject to 



and 



subject to 



dt dx \ dx 
|^(0,t) = , p(s-(t),t)= Pl (a), (8) 



|£(M)=0 , p(s+(t),t)=p 2 (a). (10) 



The evolution of s(t) is determined by the Rankine-Hugoniot jump condition 

ds _ ( J+ - J-) 



(11) 



dt (p2 - pi) 

where the fluxes are given by J ± = —D(p(s ± {t),t))d x p(s ± (t) : t), this equation being obtained by 
differentiating the statement of conservation of mass 

(*) /•! 

p(x,t)dx + p(x,t)dx = M, (12) 

Js(t) 

and using [[71)-(jT0]). 

The initial datum for p, satisfying the boundary and jump conditions, is chosen to be smooth 
away from the initial discontinuity, with p < p b in the low-density phase and p > p" in the high- 
density phase. By construction, D(p) is then initially positive on each phase, and, by virtue of the 
Maximum Principle, we can reasonably expect the Stefan problem to be well posed. In the sequel, 
we will refer to the coupled system ([7l)- ([TT|) as SP\\ analogous problems SP n , with n jumps, will be 
treated in Section 3. 

Finally, note that we allow for the possibility that s(t) may hit the domain boundary at x = or 
1 in finite time. In this case, provided the gradient remains bounded as the boundary is approached, 
the solution can be continued via the ordinary Neumann problem for (TTJ), which we will call NP. As 
a convention, a solution continued in this way will still be referred to globally as a solution of SP%. 



2.2 Steady states and their stability 

If the total mass M satisfies p\ < M < P 2 then there is precisely one discontinuous steady-state 
solution of SPi , given by the step function 

p*(x) 

s* 



pi : xe[0,s*) 
p 2 ■■ x G (s*, 1] 

(P2 - M) 

(P2 - Pi) ' 



(13) 
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while if M < p\ or M > pi then no such discontinuous solution exists. 
Moreover, SP% clearly has the uniform steady-state solution 

p{x) = M, < x < 1 (14) 

if and only if M < p b or M > pK 

For the two possible kinds of steady state we have the following stability results: 

Theorem 2.1 Suppose we have a global smooth solution pair, (p(x,t),s(t)) } for SP\, with initial 
data (po(x), sq), such that < sq < 1. Then 

(i) if $ < M < p", p(x,t) converges exponentially to p* in L 2 -norm, and s(t) converges expo- 
nentially to s* , as t — > oo, 

(ii) if p\ < M < p b (resp. ffi < M < p 2 J> an d Po ^ M (po > M) in the low (high) phase, 
the solution is attracted towards (p*,s*) as in (i), 

(Hi) if M < pi (resp. M > pi), then s(t) hits x = (x = 1) in finite time, and the continued 
solution converges to the uniform steady state p = M, exponentially in L 2 . 

(iv) if M — pi (resp. M — pi), then s — + 1 (0) as t ST, some T, with the possibility that 
T = oo, and p(x,t) approaches p\ (p%) exponentially in L 2 in the low (high) phase as t S T. If T 
is finite, the continued solution of NP converges to the uniform steady state p = M , exponentially 
in L 2 . 

Proof. To prove (i), first note that, by (fT2|) and the fact that p cannot enter I a = (p b ,p^), 

< s min < s(t) < s rnax < 1 ,Vt, (15) 

where 

Smin(M) = {p*-M)/ P \ (16) 
Smax(M) = (l-M)/(l-p b ). (17) 

Next note that, on the subinterval (0, s(t)), we have 

Hence, multiplying through by (p ~ pi), integrating by parts, and using the mixed boundary condi- 
tions gives 

i a /- s (*) /■«(*) 

' ' (p- Pl ) 2 dx = - D(p)(d x (p~Pi)) 2 dx, (19) 



2dt 

which implies 

~-t;\\p- Pi\\h{o,s) < -4dx{P~ ^i)lli=(o, s ). ( 2 °) 



(i 



2dV 

where e = ini 0<x<s{a) {D(p(x, 0))} > 0. 

Since p = p\ at x = s~(t), we also have the Poincare inequality 



\\p - Pih^o.s) < s\\d x (p - pi)||i2 (0)S ), (21) 
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and hence (f2T))) implies 



|]p-Pi|lz,=(o, s )(*) < \\P- Pi\\h(o,s)W e 



(22) 



An analogous inequality on the subintcrval (s(t), 1) is obtained in exactly the same way. 
Next, from (p~2|) we have 




/ (p(x,t) - px) dx + pis(i) + / (p{x,t) - p 2 ) dx + p 2 (l-s(t)) 



= M, 



(23) 




and hence, rearranging and using the L 2 -decay just shown, 




-c 2 t 



(24) 




(25) 



as required. 



To prove (ii), note, for example, that if pi < M < p b and po(x) < M for < x < sq then p = M 
is a supersolution in the low phase. From this, it follows once again that s(t) remains bounded 
away from and 1, by conservation of mass. Convergence to (/?*, s*) is proved as before. The same 
argument goes through for < M < p 2 and po{x) > M, so < x < 1. 

For (hi), if M < p\ (resp. M > p 2 ) then the only candidate for a steady-state solution is given 
by p = M and s = 1 (s = Q), and, by L 2 -decay in the low (high) phase and conservation of mass, 
s(t) is forced to hit the boundary in finite time; _L 2 -convergence for the subsequent NP is proved 
via the usual energy estimate 

Finally, if M = pi (resp. p 2 ), then the only possible steady state is p = p\ (p 2 ) for < x < 1. 
Clearly, s(t) remains bounded away from (1) and therefore p{-,t) decays to p\ {p 2 ), exponentially 
in L 2 on (0, s) ((s, 1)), as we have already seen. The L 2 -decay implies that s — > 1 (0) as t increases, 
but the convergence may take infinitely long in this exceptional case; if the convergence takes place 
in finite time then the subsequent NP is as in (iii). This proves part (iv) □ 

We conclude this discussion by noting that, for the cases not covered by Theorem 2.1, namely 

1. pi < M < p , po{x) > M for some x < so, 

2. < M < p 2 , po(x) > M for some x > s , 

we have not been able to find a clean analytical criterion for determining which of the two possible 
steady states will be approached at large times, for given initial data. However, in the next subsection 
we show that the discontinuous steady state (|13p is always linearly stable in a certain sense, whenever 
it exists, and, moreover, we report on numerical simulations which suggest that s(t) can hit the 
domain boundary in finite time, provided condition 1 or 2 (directly above) is satisifed, and the 
initial density profile is sufficiently far from (Ti"3"|) ; in other words, SPi appears to be bistable for 
some values of the total mass. 
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2.3 Existence and uniqueness for smooth data 
2.3.1 Background 

We can make a connection between our moving-boundary problem SPi and the traditional two- 
phase Stefan problem for the melting of ice in water by identifying p as the specific latent heat, and 
a := K(p) as the temperature. In the traditional formulation, our low-density phase is thought of 
as the solid, and our high-density phase as the liquid phase. Given this, one imposes the physically 
reasonable condition 

Po < Pi m the solid (or low-density) phase, po > p2 in the liquid (or high-density) phase, (26) 

and then one can write down a favourable weak formulation of the problem, which can be shown to 
have a unique solution for bounded data (see Appendix C). Furthermore, this weak solution turns 
out to be a global- in-time smooth solution of the original problem (see, e.g., [5]). 
For SP\, however, the restriction on the initial data is just 

Pa < p b in the low-density phase, po > p" in the high-density phase, (27) 

which is weaker than (|26[) . In particular, this entails that K(p) can no longer be assumed mono- 
tonically increasing, and, consequently, that there is no nice weak formulation of the problem. We 
note that, in the literature, uniqueness of solutions, given (|26j) . is usually proved using the weak 
formulation [H[5]. Moreover, in the basic proof of local existence of smooth solutions as presented in 
[5], for example, it is not so easy to discern whether the restriction (|2l)|) can be relaxed to (|27[) . For 
these reasons, we now present a concise, reasonably self-contained, classical existence-and-uniqucncss 
theory for SPi, with data subject merely to (|27|) . which employs the standard Holder and Sobolev 
estimates of linear parabolic theory. 



2.3.2 A local-existence theorem 

In order to solve the moving-boundary problem S Pi, we employ, in each phase, the standard trick 
of rescaling the spatial variable (see, for example, [7]), such that ([7]) and (fT0|) become a pair of 
fixed-boundary problems, coupled by the (rescaled) Rankine-Hugoniot condition. 
Specifically, for < x < s(t) we introduce the coordinate transformation 

(x, t) -» (x, t) : x = -, t = t, (28) 

while for s(t) < x < 1 we take 

(x, t) ->• (x ,t) : x=^—^-, t = t. (29) 
1 — s 

In terms of these new coordinates, and writing v{x,t) — p{x,t) for < x < s(t), w(x,t) — p(x,t) 
for s(t) < x < 1, equations (O and (fT0|) become, upon dropping hats, 

dv Id / )® V \ + * ® V (30) 
dt s 2 dx \ dx J s dx ' 

for x £ (0, 1), subject to v x (0,t) = 0, v(l,t) = pi, and, respectively, 

dw 1 d f dw\ . (l-x) dw 

It = (T^5 to r H to ) + s (T^y to ' (31) 
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for x £ (0, 1), subject to w(0, t) — p 2 , ^(l, t) = 0, while the jump condition assumes the form 



ds = _ / D(p 2 )w x {0,t) _ D( Pl )v x (l,t) 
dt \ (1 - s) ' s 

and we also have the rescaled conservation-of-mass equation 



(P2 - Pi) 



Al- 



io w dx 



f Q v dx — f w dx 



(32) 



(33) 



Finally, for the application of classical parabolic theory, the initial data and the Dirichlet con- 
dition at Sq :— s(0) must satisfy a certain first-order compatibility condition, which is obtained by 
setting v t and w t equal to zero in (|30p and (|31[) . and substituting the initial data, evaluated at the 
phase boundary, into the right-hand sides of (|3"0|) , (|3"Tj) . and (|3"2"|) . 

The most compact way of writing this compatibility condition is to go back to the original 
coordinates and introduce the dependent (temperature) variable a = K(p), in terms of which the 
relevant equations are simply 



(34) 



where, for example, = <r x (s ± (0), 0). 
We are able to prove the following: 



Theorem 2.2 Given initial data (vo,wq, sq) , such that vo,wq G C' +2 ((0, 1)), I > 0, < sq < 1, 
< Vq < p^ , ffi < Wo < 1, and such that the first-order compatibility condition is satisfied, the 
system of equations i30\) , i31\), \33\) (and hence problem SP\) has a unique classical solution on 
some small time interval [0, T] . 

Proof. The proof proceeds via an iterative scheme of successive approximations, and to get the 
required strong-convergence properties, we work in the setting of 'parabolic' Holder spaces. 

First of all, we set Qt = [0, 1] x [0, T], and, as in [4], [5], for example, let H M/2 (Q T ), I > non- 
integer, denote the Banach space of functions u with continuous derivatives D r t D s x u, for 2r + s < I, 
equipped with the norm 



i(0 



['] 



= £ \\d;d' x u\\ L o B(Qt) + J2 {DlDiut;^ 

2r+s=[l] 



(DlDlu)l Q l 



s v(2-2r-s)/2 



2r+s=0 

where, for < [3 < 1, 



0<l-2r-s<2 



(35) 



sup { \v(x, t) — v(x', t)\ \x — x'\ l3 } , 

(x,t),(x',t)EQ T 

sup {\v(x,t) - v(x,t')\\t - i'l^} . 



(36) 
(37) 



Now we describe our iterative scheme for obtaining new approximate solutions of (f30|) . (j3T|) . ([32]) 
from old ones. 
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Suppose we have smooth ith iterates (v l ,w l ), with positive diffusivities on Qt, which satisfy the 
initial and boundary conditions, and suppose that 

m^{\v%l 1 \\w% + T 1) )<C, (38) 

for some I £ (0, 1). 

Next, determine the approximant s'(i) by solving the nonlinear ODE 

ds* ( D(p 2 )wj(0,t) D( Pl )vUl,t) \ , x 

^ = -{ 7 , (39) 

which, for the given datum sq, has a unique smooth solution on [0, T], provided T — T{C) is chosen 
small enough (Picard's Theorem). Clearly, the ordinary Holder norm of s l is bounded according to 

II s \\C l ' 2 ([0,T]) < F(C), (40) 

for some function F. 

The functions v l , to' , s l are now inserted into the right-hand sides of the linearised field equations 
dv i+1 1 d ^dv i+1 \ s*<9u i+1 

, (41) 



dt (s 1 ) 2 dx \ dx J s l dx ' 

for x £ (0,1), subject to v t+1 (- 7 0) = v , (v i+1 ) x (0,t) = 0, v i+1 (l,t) = p 1} and 

dw t+1 1 d / , i,dw i+1 \ .,(l-x)dw l+1 



dt (1 — s 1 ) 2 dx \ dx J (1 — s l ) dx ' 

for x £ (0,1), subject to w l+1 (-,0) = w , w l+1 {0, t) = p 2 , {w t+1 ) x {l, t) = 0. 

These equations are parabolic, by the choice of (v % , xif) , and therefore have unique smooth solu- 
tions v l+l , w l+l on Qt, by Lemma A.l (see Appendix A), which is a convenient summary of relevant 
results from classical parabolic theory [4]; the updated approximants also have corresponding posi- 
tive diffusivities, by the Maximum Principle. 

The coefficients of (j4l"j) . (j42|) . when expanded into the standard form (|85|). are such terms as 
xs l I V, D(v l ) I '(s 1 ) 2 , D' (v l )v l x / {s 1 ) 2 , and analogues for w\ and are therefore dominated in H 1 ' 1 / 2 (Qt) 
by and Thus, by Lemma A.l, 

K +1 |£ 2) <F(C)(K.,0)|^+ Pl ), (43) 



and 



l^ +1 |^ +2) < 



F(C)(K-,0)|[; +2) -t-p 2 ), (44) 

where F is some positive function, il = (0, 1), and | • \q +2 ^ is the ordinary Holder norm. 
Next, applying Lemma B.l to (j43|) and l[44|) results in 

-if(n)|| Uo || c2 < *T(n)T*F(C) (|« |g +2) + P1 ) , (45) 

and 

-^(flJIKIIo- < X(«)T 5 F(C) (H)|£ +2) +P2) , (46) 

where S = min {4 > f (1 — 0} 
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Thus, if we choose C > K(d) max(||v(-, 0) ||c 2 , 0)||<7 2 )> and T{C) is taken sufficiently small, 
then 

max(|^ 1 |g 1 )>^|^)<C7. (47) 

Iteratively, we therefore have that ([25]), (g0|), (TJ3|), and (gj) hold uniformly for all i. 
Next, taking the difference of (|4"Tj) and the corresponding equation for v % gives an equation of 
the form 

J^+i - „<) = Di (x, t)^(v l+1 ~ v l ) + - v*) + bi(x, t)-j^{v< ~ v*- 1 ) 

+ c l (x,t)(t; I -w l - 1 ) + d l (x,t)(s I -s l - 1 ) + e l (2;,i)(s l -s 1 - 1 ), (48) 

11/2 

where Di > e for all i, some e > 0, and all the coefficients are bounded in Hq^ , uniformly in i, by 
the results just obtained. 

For the difference s l — s 1 " 1 , we get, from (|32|) . 



|( fl < _ s i-i) = p.(i)( fl < _ + _ „i-i)(i, t ) + r . _ wr i)( , t ), (49) 

where pi, qt, Ti are bounded in C l l 2 , uniformly in i, and from this it is easy to deduce that 

\ A i _ < d (V - + K - ta*- 1 !^) . (50) 

Thus, applying ([50]) and Lemma A.l to (|4*5|) . we get 

K+ 1 - v% + V < C X (y - v*- 1 ]^ + K - to*- 1 !^) . (51) 
Of course, there is also an analogous estimate for w l+1 — w l , which, together with (fSTj) . implies 



i« t+i - «t: 2) + k +1 - <i 2) < & (V - + k - ■ (52) 



Finally, an application of Lemma B.l to the right-hand side of this inequality gives, since all 
iterates have the same initial data, 



< C,T S (y-v*-^ + k-^" 1 ^) . (53) 



(Note that C\ stands for various constants in the above). 

If T is chosen so small that C\T < 1, it is easy to see that v l and w l are Cauchy sequences with 
respect to | • |g^ 2 - ) -norm. Moreover, it follows from (|50"|) that s l is Cauchy in C 1+ ^-norm. By the 
strong convergence of these sequences, their limits v, w, s satisfy (|30|) . (|3"Tj) . (|32|) pointwise, and thus 
constitute a classical solution of SP\ . 

Turning to the question of uniqueness, suppose we have two smooth solution triples (v,w, s), 
(v,w,s), satisfying (|30[) . (|3ip . (|32[) . and having the same initial data. Subtracting the equation 
satisfied by v from that satisfied by v then gives an equation of the form 

d t (v - v) = D{x,t)dl{v -v) + a(x,t)d x (v -v) + b(x,t)(v - v)+c(x,t)(s - s) + d(x,t)(s - I), (54) 
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where the coefficients are smooth, and D(x,t) > e, some e > 0. An analogous equation is obtained 
for w — w in the same way. 

Thus, by Lemma A.l, v — v can be bounded by the inhomogeneity in (|54p . according to 

\v-v\ [ q t 2) < c(\c(xM*-¥$ T + \d(x,t)(i-¥$ T ] 

< C (\s-s\g + \s-S\g) 

< C (\v x -v x \q t + \w x - w x \^ 

< c(\v-v\^ + \w-w\^), (55) 

for various constants C, where we used the Rankine-Hugoniot condition to get the third line. 
Adding the analogous inequality for w — w therefore results in 

\v - + \w- <c(\v- v\« + V + {w _ ^|gD) ( (56) 

which, with the help of Lemma B.l, implies 

I- ~ -lg: 2) + I- - -l£ 2) < CT* (\v - + \w - w\ ( ^) . (57) 

Choosing T small enough gives, finally, v{x, t) = v(x,t), w(x,t) = w(x,t), as required □ 

2.3.3 Global existence subject to a sign condition on the solution gradient at the 
discontinuity 

We will prove a global-existence theorem for SPi by showing that, for a local classical solution 
(p(x, t), s(t)), the gradient, p x {x,t) 1 and hence also s(t), are a priori bounded, provided that the 
one-sided limits p x (s ± (t) , t) satisfy a sign condition. The proof is rather different from, and somewhat 
shorter than, that presented in [5]- 

First, with the change of dependent variable a = K{p), which is smooth and invertible in each 
phase, the governing PDE becomes 

a t = D(a)a xx , (58) 
where D(a) = D(K^ 1 (a)) > 0, the boundary conditions become 

o-x(O) = tr B (l) = 0, a(s-) = K{ Pl ) := ax, a{s+) = K{p 2 ) := a 2 , 
and the Rankine-Hugoniot condition takes the form 

ds _ (<t x (s+) - <r x (s-)) 



dt (P2-P1) 
Next, note that integration by parts gives, for any smooth a, and m = 1, 2, . . ., 

T (t) d t {a x f m dx = - f ^ 2m(m - l)a 2 x m - 2 a xx a t dx + [Zmof-Vjf 1 
Jo Jo 



(59) 



(60) 



Assuming now that (a, s) solves SP\, we have o~(s (t),t) = <r\, and hence, by differentiation, 
a x s + at — at x — s _ , which, with the aid of ([SO]) , leads to 
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\2m • 

I X — S~ 



^ll(^ril^(0, s ) = - f 2m(2m-l)(<T x ) 2m - 2 <j xx <j t dx + (1 - 2m)(a x f m s 

< (l~2m) S (t)(<T x ) 2m \ x=s -, (61) 

where we used ([55)) to discard the integral term. 
By a similar calculation, there also follows 

|llK) m H W) < -(! - 2m)^)M 2 "W- (62) 
If we regard o x as a function in L 2 (0, 1), then adding (IFl]) and (f6"2")> . and using (|55|) . results in 

4l(^ril5 < (1 - 2 m ) iaAs p- a ^ ) \ (a x f m (s+) (* x ) 2m (s-)). (63) 
di (/9 2 - Pi) 

Now, if 2m = 2", n = 1,2,... then, inductively, there is a positive multinomial F n (a, b) such 
that 



(a -6) (a 2 " -6 2 ") = F n (a, 6)(a + 6), Va, 6. (64) 

Using (|64|) in the right-hand side of (|63p with a = o- 2; (s + ),6 = a^s - ) shows that ||(ca;) m ||2 is 
decreasing for 2m = 2™, provided 

(r x (s+) + a x (s~) > Vi. (65) 

This condition can be guaranteed, for example, by choosing data (po,so) f° r 'S'-fi satisfying the 
traditional condition (|26p . and by appealing to the Maximum Principle in each phase. 
Thus, assuming (|65[) holds, we have 

< Ikxll^-(o) < Ikxlloo(o), (66) 

for n — 1, 2, . . ., and therefore, by Theorem 2.8 of pQ, 



IkxIlooW < IkxlU(O). (67) 

In 1-d, such an a priori bound on |j Vcr||oo, which of course also gives a pointwise bound on the 
corresponding p x and on s(t), is actually enough to continue the local solution of SPi obtained in 
Sect. 2. 3, by standard theory. Indeed, if v and w are as in ([30 |) -(f3"T j) . then combining the estimates 
of Appendix A gives, for 3 < q < 4, and a constant C which is controlled by Va, 

Mq'^ <C(\\v\\ 

c 2 ((o,i)) + Pi) i (68) 

and an analogous estimate for w. 

Thus, for some (3 > 0, the H@ +1 ' a ^ +1 ^-norm of v and w is controlled by Va and the C 2 -norm 
of the initial data. Given this, Lemma A.l implies in turn that the H^ +2 ' 3 (^ +2 )-norm of v and w 
is controlled by Va and the C^ +2 -norm of the initial data. Thus, the local classical solution can 
always be extended onto to a larger time interval, provided, of course, s(t) stays away from the 
domain boundary. 

We have already seen (Theorem 2.1) that if the total mass M lies outside the interval (p , p"), 
it cannot be ruled out that s(t) will hit the domain boundary in finite time; should this occur, the 
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solution can be continued via the Neumann problem on = (0, 1), as mentioned above. Note that, 
on a disappearing phase, the solution (in terms of the original variable p) merely approaches p\ or p2, 
as appropriate, by the boundedness of the gradient. Since we have already extended our definition 
of SP\ to cover such eventualities, we have therefore proved 

Theorem 2.3 Given phase-wise smooth initial data (po(x), sq) satisfying i26\) and the first-order 
compatibility condition, problem SPi has a unique smooth, global-in-time solution (p(x,t), s(t)), such 
that the corresponding Hct^Hoo is monotonically decreasing for allt. 

In general, it is not clear whether a corresponding global-existence result can be obtained in the 
situation where (|2l))) is not satisfied. An exception to this is the special case of 1-phase problems, to 
which we now turn our attention. 



2.3.4 Global existence for 1-phase problems; Mass Lagrange coordinates 

If the initial density is constant for either x < sq or x > sq, then we refer to the corresponding 
evolution problem SP\ as a one-phase problem. As a consequence of the Neumann condition at the 
domain boundary, it then turns out that there is a change of the spatial variable, different from a 
simple rescaling, which transforms the one-phase SP± into a regular quasilinear parabolic problem 
on a fixed domain. This allows us to prove a global-existence theorem, using standard parabolic 
theory, regardless of the direction of the gradient at the phase boundary. 

Without loss of generality, let us assume that p is constant (= p<z) in the high-density phase, and 
variable in the low-density phase, and let us again make the change of dependent variable a = K(p) 
on < x < s(t). Since we are assuming p < p b in the low-density phase, this change of variable is 
invertible there, with inverse denoted by p = b(a) := K (a). 

In terms of a, the governing equation {TJ) becomes 

-¥ = ^ (69) 

subject to o~ x = at x = 0, and a = K{p\) at x = s(t), while the Rankine-Hugoniot condition takes 
the form 

ds da , . ., , , . 

* = (70) 

Next, we transform by introducing so-called Mass Lagrange coordinates [5], (r,y), which 
are defined by 

s(t) 

r = t, y= J [p 2 -b{a{x,t))\ dx. (71) 

x 

Since, by construction, ^ = — (p2 — b(o-)) < 0, this coordinate change is good, and, introducing 
v(y,r) — a(x,t), ([55]) now takes the form 



b'(v) dv d 2 v 



(p 2 - b(v)) 2 8t dy 2 ' 
which is a regular quasilinear parabolic equation for v. 



(72) 
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Clearly, the moving boundary x = s(t) gets mapped to y — 0, and, as a consequence of conser- 
vation of mass, the domain boundary x = gets mapped to 

S (0) 

2/o := J [p 2 — b{<r{x,0))) dx — const.. (73) 
o 

Thus, ([7^)1 is to be solved on a fixed spatial domain subject to the boundary conditions || = at 
V = Vo, and v = K(p 1 ) at y = 0. 

It is well known that this problem has a unique, global classical solution v(y, r), given H 1 initial 
data compatible with the Dirichlct condition at y = [5J[1]. The corresponding jump location s(i), 
which can be reconstructed from v(y, r) by integrating, for each t = r, the equation 

^=b(v(y,r)) (74) 

from the point (x = 0, y — yo) until y hits zero, could conceivably hit x = 1 in finite time (see the 
numerics in the next subsection). If this occurs, the solution of (|72|) should be stopped, and then 
continued for all time via NP for UJ on CI = (0, 1). 
We thus arrive at 

Theorem 2.4 Given smooth, one-phase initial data which satisfies the first-order compatibility con- 
dition, problem SPi has a unique smooth, global-in-time solution (p(x,t), s(t)). 



2.3.5 Linear stability near a discontinuous steady state 

The use of rescaled coordinates, as introduced early on in this subsection, allows us to investi- 
gate linear stability of the unique discontinuous steady state (fTB"]) of SPi, which exists as long as 
Pi < M < p 2 . 

Linearising (|3"0|) and (|3ip around v — p\ and w = p 2l respectively, and using hats to denote 
differentials, we get the pair of heat equations 

«i = T^y D (Pi)vxx, (75) 

w t = ^ — s ^ 2 D(p 2 )w xx , (76) 

which are to be solved subject to v(l) = 0, ^(O) = 0, w(0) = 0, w x (l) = 0, while linearising ([33]) 
around s = s* gives 

s - (Pi - P2Y 2 -Pi)f w dx + (p 2 - M) (v - w) dxj. (77) 

For a solution (v, w, s) of the linearisation, it is therefore clear that, say, ||u||oo — * and ||u)||oo — * 
as t — > 00, and also that s — > as t — > 00. In this sense, then, the unique discontinuous steady state 
of SPi is always linearly stable, whenever it exists. 



2.4 Numerical simulations 

As well as facilitating mathematical analysis, the rescaled, fixed-boundary representation of SP±, 
(|3H)) . PT]) . (|3"!?|) . also comes in useful for numerical simulations. 
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(a) 



(b) 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Figure 2: Solution of SPi with a — 0.85 and initial data satisfying (|26|) . Density profile shown at 
(a) t=0, (b) t=0.0125, (c) t=0.1289, (d) t=0.456. 

Specifically, we use a method-of-lines approach in which the diffusion terms in (|30p and (|3ip are 
discretised using the random- walk model ([3]) , while the advection terms are discretiscd by means of 
a standard, explicit upwinding scheme. A simulation in which the global-existence criterion (|26[) is 
satisfied is shown in Figure 2, and one in which it is violated is shown in Figure 3. In each case, 
the solution approaches the appropriate discontinuous steady state at large times. Several other 
simulations have been carried out in the case where (|26p is violated, and no singularities have been 
observed to develop. 

Turning to the question of bistability for cetain values of the mass M, we next show in Figure 
4 a close-up of a solution with M = 0.3184 and a = 0.85, for which the phase boundary hits the 
domain boundary in finite time. This should be compared with the simulation shown in Figure 5, in 
which the solution has the same mass and adhesion coefficient, but this time evolves to the unique 
discontinuous steady state as t ~ -» oo. 

It is perhaps worth noting that such bistability cannot occur if the traditional restriction (|26l) is 
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Figure 3: Solution of SPi with a = 0.85 and initial data violating (j26|) . Density profile shown at 
(a) t=0, (b) t=0.0125, (c) t=0.1289, (d) t=0.456. 
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imposed on the initial data, by Theorem 2.1. 



3 Multi-phase problems 

We now consider the general case of our Stefan-problem set-up, in which the initial density profile 
jumps n times across the the unstable region I a . This initial/boundary- value problem will be denoted 
by SP n . 



3.1 Formulation 

The problem SP n consists of looking for a piecewise-smooth p which satisfies the diffusion equation 
(fl]) away from n discontinuities Si(t), with < Sj < s^+i < 1, such that p always jumps between pi 
and p2 at the Si, which evolve according to the Rankine-Hugoniot condition 



dsi {Jt-Jl 



dt [pi 



(78) 



where jf = ~D{p(sf(t),t))d x p{sf(t),t), and [ Pl ] = (p(s+ ',*) - p( S r ,t)) (= ±(p 2 - pi)), and the 
Neumann condition is again imposed at the domain boundary x = 0, 1; a schematic for SP n is given 
in Figure 6. 

We also allow for the possibility that s\ or s n may hit the domain boundary in finite time, or 
that neighbouring discontinuities could collide, leading to the annihilation of a phase. Should any 
of these events occur, the solution can be continued via SP n -\ or SP n -2, as appropriate, and so on. 
A solution continued in this way will still be referred to globally as a solution of SP n . 



3.2 Steady states and their stability 

The steady-state picture for SP ni given n > 2, is more complicated than that for SP\. 

First of all, for a given mass M satisfying p\ < M < p2, there is a continuum of two- valued, 
n-jump steady states, each of which is given by a choice of the Si which merely has to be compatible 
with M. Moreover, discontinuous steady states with fewer than n jumps can also be considered 
permissible - these could be approached dynamically by (multiple) coalescence events, and/or by 
(successive) collisions of phase boundaries with the domain boundary. Each steady state with more 
than one discontinuity is expected to be merely neutrally stable, since the total mass is invariant 
under small translations of an internal phase (i.e., one which does not touch the domain boundary). 

If the stronger condition p b < M < p" holds, then there is no possible continuous steady state, 
and thus, for a global solution of SP n , at least one discontinuity must remain as t — > oo. If, instead, 
pi < M < p b or < M < p2, then the uniform steady state exists alongside the discontinuous 
family already discussed. In this regard, note that Figures 4 and 5 could be considered as simulations 
of one half of a reflection-symmetric SP2, with pi < M < p b , in which the central high-density phase 
is either annihilated in finite time (Figure 4) or preserved as t — > 00 (Figure 5), depending on the 
proximity of the initial data to the (unique) discontinuous, symmetric steady state. 

In general, it is difficult to say anything analytical about the stability of steady states of SP n , 
although in each phase the solution will still decay to p\ or P2, as appropriate, exponentially in L 2 , 
by essentially the same calculation as in Section 2.2, as long as the relevant phase boundaries stay 
away from the domain boundary. 

Finally, note that if M > P2 or M < pi, then there is only one possible steady state, namely the 
uniform one, and, by L 2 -decay, phase boundaries for a global solution must disappear in finite time 
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(a) (b) 




0.7 0.75 0.8 0.85 0.9 0.95 1 0.7 0.75 0.8 0.85 0.9 0.95 1 




0.7 0.75 8 0.85 0.9 0.95 1 0.7 0.75 0.8 0.85 9 0.95 1 



Figure 4: Close-up of a solution of SP\ with M = 0.3184 and a = 0.85, such that s(t) hits x = 1 in 
finite time. Density profile shown at (a) t=0, (b) t=0.001, (c) t=0.0051, (d) t=0.012. The density 
to the left of x = 0.7 is essentially constant throughout the simulation. 
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Figure 5: Solution of SP\ with M = 0.3184 and a — 0.85 which evolves to a discontinuous steady 
state as t -> oo. Density profile shown at (a) t=0, (b) t=0.0373, (c) t=0.6264, (d) t=0.7479. 
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8,(1) S 2 (t) 8,(1) S i+1 (t) S„(t) 1 X 

Figure 6: Example of a density profile for the multi-jump Stefan problem SP n - 



(or possibly infinite time in the exceptional cases M = p\, P2) via coalescence events, or by merging 
with the domain boundary. For the subsequent NP, exponential i 2 -convergence to the uniform 
steady state p = M follows as before. 

3.3 A local-existence-and-uniqueness theorem for smooth data 

The multi-phase problem, SP n , is solved, locally in time, in the same way as SPi; in each phase 
the spatial variable is rescaled in order to fix the moving boundary (or boundaries), and the same 
estimates go through as before. 

To be explicit, on the interval [sj, Sj+i] we make the coordinate transformation 



(x, t) -> (x, t) : x = - -, t = t, 

[Si+i - s l ) 



(79) 
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such that, writing Vi{x,t) — p(x,t), the governing PDE ([T]) becomes, upon dropping hats 

~dF = {Si+i ~ Si) d-x\ D{Vi) ^) + (s i+1 - Si ) to' (80) 

on [0, 1] x [0, T], for i = 0, 1, 2, ... , n, where we adopt the conventions so(t) = and s n+ i(t) = 1. 

On an internal phase, this equation is to be solved subject to the appropriate Dirichlet conditions 
at x = 0, 1, while on an extremal phase the mixed Neumann/Dirichlet conditions are used, as in the 
two-phase case (see Figure 6). 

The rescaled family of Rankinc-Hugoniot conditions takes the form 

dsi = _ / D(p2)(v i+1 ) x (p,t) _ D{pi)(vj) x (l,t) 
dt \ (sj+i - Si) (si - Si-i) 

for i — 1, 2, . . . , n. 

Equations (|80|) . (|8T|l are solved by the same kind of iteration employed for SPi; the required 
Holder estimates are obtained via linear parabolic theory and Picard's Theorem for systems of ODEs. 
The upshot of all this is 

Theorem 3.1 Given initial data W®, s°), i = 1, 2, . . . , n, for which the and w® belong to the 
Holder space C l+2 , I > 0, and such that the first-order compatibility conditions are satisfied, the 
system of equations t80\) , &81\) , and hence problem SP n , has a unique classical solution on some 
small time interval [0, T]. 

3.4 Continuation of the local solution, subject to a sign condition on the 
solution gradient at discontinuities 

With the same notation as in Sect. 2.4, and by a similar calculation, it is straightforward to see that 
the gradient of a local solution of SP n satisfies, for m = 1, 2, . . ., the a priori estimate 

|llK) m lli 2 ( ( o,i)) < Ofzf) E si s n ^ M«< + ) - °*(»T)) (°T l (4) - <£"(0) • (82) 

V P2 Pi / ^ 

Thus, Halloo, and hence also the ||px||z oo (( s «,si + i)) an d Sj, are bounded for all time, provided 

sign[p]i (a x {s+) + a x (s7)) > 0, Vi. (83) 
This condition holds if, for example, 

Pa < Pi in low-density phases, po > P2 in high-density phases, (84) 

which is the multi-phase analogue of (|2U|) . 

Inequalities ( 84) therefore guarantee global existence of the corresponding solution to SP„ . mod- 
ulo coalescence events, and the possibility that an extremal discontinuity might hit the boundary in 
finite time. 



)(P2-Pi)-\ (81) 
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Sift) s i+1 (t) 



Figure 7: Two coalescing high-density phases; note that the resulting single high-density phase has 
a corner at the point and time of coalescence. 
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3.5 Coalescence events, and continuation thereafter 



In order to continue a solution of SP n after the coalescence of two phases, it is necessary to extend 
our existence theory (which has thus far required the initial data to lie in C l+2 ) to the case where 
the initial density profile may have a 'corner' in one of the phases. This can be seen by considering 
the situation illustrated in Figure 7, in which two high-density phases coalesce, thus annihilating a 
low-density phase. 

If we make the assumption that condition (184|) holds, then, in fact, the required extension for 
such phase- wise H 1 data follows from Theorem 14, Sect.V of [5], and hence, putting together the 
results of this section, we therefore arrive at 

Theorem 3.2 Given n initial jump locations (si)o, and a phase-wise smooth initial density profile 
Pq satisfying the first-order compatibility conditions at phase boundaries, along with the gradient 
condition £84]) , the multi-phase Stefan problem SP n has a unique, global classical solution, on the 
understanding that phases may in time be annihilated via coalescence events, or that an extremal dis- 
continuity may hit the domain boundary infinite time. Moreover, ||(T x ||oo is monotonically decreasing 
for all time. 

We end by noting that it is not clear whether even a local existence theorem for SP n can be 
proved for merely -ff *-data if (|84|) is not satisfied; one can, for example, proceed by approximating 
H 1 data with smooth data, but it seems that the a priori bound on a x implied by ([84]) is required 
to get the necessary convergence. 

4 Concluding remarks 

The results of this paper can be interpreted as saying that cell-cell adhesion is enough to stabilise 
sharp-edged aggregations of diffusing cells if there is enough mass present in the biological domain, 
and that diffusion homogenises the cell density in each given high- or low-density phase. For low 
masses, however, narrow high-density regions can become annihilated even if the adhesion is very 
strong - that is to say, diffusion in the low-density regions wins, leading to a globally uniform cell 
density at large times. In the case of intermediate masses, bistability becomes possible - either 
adhesion or diffusion can win out, depending on the profile of the initial data. 

Next, while being motivated by the adhesion-diffusion equations (fTJ) and ©, it should be clear 
that the analysis carried out in this paper does not require that the diffusivity D(p) have the special 
form (J2]); indeed, all arguments go through for any equation of the form p t — K{p) xx , provided the 
C 2 function K is increasing outside an unstable interval of p values, and provided [K(p)]i — at 
jump locations Si. 

Finally, one rather obvious biologically-relevant extension of the work described here would be 
to carry out a similar analysis with a chemotactic term factored into the right-hand side of ([T]). In 
that situation, one could imagine beginning with a low-density initial datum, evolving the solution 
until, through chemotactic aggregation, p hits the unstable region I a at some point x c , and then 
continuing the solution via a Stefan problem with an initial spike at x c which jumps from p\ to pi. 
The transition from well-posed Neumann problem to Stefan problem is somewhat singular in that 
case, and is the subject of ongoing analytical investigation. 
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Appendix A. Linear parabolic theory 

First of all, we have 

Lemma A.l The equation 

du ii du 

— = D{x, t)-g^+ a(x, t)— + b(x, t)u + f{x, t), (85) 

with initial data uq G C' +2 ([0, 1]), < I < 1, subject to the boundary conditions u x (0,t) = 0, 
u(l,t) = u\, and the first-order compatibility condition 

D(x, 0)^? + a(x, 0)^ + b(x, 0)u + f(x, 0) = 0, (86) 
has a unique solution u G H l+2 ^ l+2 ^ 2 (Qx) on Qt — [0, 1] x [0, T], which satisfies the estimate 

(l/lg T + K||c«+"([o,i]) + l«i|), (87) 

with the constant c remaining bounded as T — > ; provided D(x,t) > e > and all the coefficients 
and the inhomogeneity in 185\) are bounded in H l ' l l 2 {VL T ). 

For the proof, see g], Ch.4, Thms. 5.2-5.4. 

Next, let H 2 ((0, 1)) denote the usual second-order L 2 -type Sobolev space on the unit interval, 
let W 1,9 ((0, 1)), q > 1, be the Sobolev space with norm 

IMk? = IMIl*((o,i)) + IWU*((o,i)), ( 88 ) 

and let W 2 ~~i' q ((0, 1)) denote the fractional-order Sobolev space with norm 

l«ll 2 -a, a = Nk, +{jj ^^g^ dxdy } , (89) 



\ut + J ] < c 




where Q = (0, 1), a = 1 - |. 

Then it is elementary to prove 

Lemma A.2 # 2 ((0, 1)) ^> W 2 ~^ q ((0, 1)), for 1< q < 4. 

Finally, let VF 2,1 (Qr) be the L 9 -type parabolic Sobolev space with norm 

h\& T = II«II,(Qt) + KIWQt) + \\u xx \\ q (Q T ) + \\utUQ T ). (90) 

Then, by the fractional-order case of \4\, Ch. IV, Thm. 9.1, we have 

Lemma A. 3 If u is the solution of the mixed Dirichlet/Neumann problem for (f8"5"j) . and if 

1 1 I |oo is a priori bounded, then we have 



(lK,o)|| 2 _2 )9 (n) + K|). (9i) 
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Finally, by gj, Ch. IV, top of p.343, we have 
Lemma A. 4 For u <E W 2,1 (Qt), and q > 3, there holds the imbedding inequality 



(92) 



Appendix B. A useful inequality 



Lemma B.l For a function u e H 2 +'<2( 2 +0 (ft T ) ; o < I < 1, such that n T = H x [0,T], SlcR, 
we have 

\u\^ < C(Q,l) (T s \u\^ 2) + ||u(., 0)11^(0)) , (93) 

where 6 = min{Z/2, (1 - l)/2}. 

Proof. The parabolic Holder norms in question are, when written out in full, 



|0+2) 



l|w||c(n T ) + IKIIc(n T ) + (u x ) l x ,n T 
\ u lt,ci T + \ u x/t,n T i 

\\u\\c(n T ) + IM[c(« r ) + \Wxx\\c(n T ) + INIc(n T ) 



+ («t)L,nr 

+ (<% T + (Ux)[%l )/2 + (Uxx)[% T . 

We proceed to estimate each of the terms appearing in |it|^ + . 
First, 



H| c(nT) <T|K|| c(nT) + || U (-,0)|| c(n) . 



(94) 



(95) 



Second, 



Third, 



||w x ||c(n T ) < swp\u x (x,t) -u x (x,0)\+swp\u x (x,0)\ 



x,t 



< T^y^u x )^ + \\u x (-,0)\\cin T) . 



(ux) l x ,n T 



\u x (x,t) - u x (x',t)\ 
sup 1 - n i 

X.x' ,t \£ X | 

\f*u xx (x",t) dx"\ 
= sup _ 

x.x'.t \*b | 

< sup|a;-a;'| 1-i [[tta;a ! ||c?(nT) 

x,x' 

= c(n,i)\\u xx \\ c(nT) 

< Csup(|u ra (x,<) - u xx (x,0)\ + \u xx (x,0)\) 

x,t 

< CT l ' 2 (u xx )% T + \\u xx {-,Q)\\c{ny 



(96) 



(97) 
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Fourth, 



, .(j+i)/2 _ HM)-mQM')I 
{)t ^ ~ S |t-tf|G+U/a 

< IK[[c(n T )T (1 - ,)/2 . (98) 



Fifth, 



j/2 _ KQM) - Mx(^,^)| 

< T^u x )^/\ (99) 
Putting these estimates together gives the desired result □ 

Appendix C. Weak formulation of the multi-phase problem 

Subject to a restriction on the initial data, a weak formulation of SP n can be written down in such 
a way that the Rankine-Hugoniot condition becomes 'hidden', thus aiding mathematical analysis. 
Essentially the same kind of formulation was written down many years ago by Ladyzenskaya et at 
4:, and, for example, their uniqueness proof goes through without change. 
For the construction that follows, we are forced to assume that 

(CI) p < pi in low-density phases, and p > P2 hi high-density ones. 
Recall that on each phase we have 

% = a*oo. ( l0 °) 

and that the Rankine-Hugoniot condition (fTTj) is satisfied at each jump. Since K(pi) = K(p2), we 
can define a somewhat flattened K(p) by 

K(p) = { „ ™ . : Pifo*} (101) 
V ' \ K{pi){= K{p 2 )) : pe[pi,p 2 }, v ; 

such that K has a piecewise-smooth, monotonically-increasing inverse, which we denote by 6; the 
functions K and b are depicted in Figure 8. 

If we introduce the new independent variable a — K(p), then clearly the equation 

^b(a) = Aa (102) 

is satisfied in each phase, by assumption (CI). 

In terms of a, the jump condition takes the simple form 

rffi = J°M) m 

dt (p2~Pl) V ' 

for i = 1, . . . , ri. 
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K(P) 



(b) 

m 



Figure 8: The functions K(p) and b(a). 

To obtain the correct weak formulation on Qt = [0, 1] x [0, T], first note that for a smooth test 
function <j)(x,t) such that <p(x,T) = 0, we have, for classical solutions of SP n , and with p := b(a), 



— / b(a)(f> dx = Si[p]i4>(si,t) + / {(j>d t b{(i) + b{a)d t 4>) dx 
Jo i= j Jo 

/ b{a)<l>{x,a) dx = 1 Ts i [p] i <l){s i ,t)dt+ / {4>d t b{a)+b{a)d t 4>) 
Jo Jo „■_-, Jo Jo 



(104) 



dxdt, (105) 



where [p]i is the leap of p at s^. 

For the right-hand side of (|102p . we have the weak form 



pT pi pT pi pT n pT n 

/ / c/)Aadxdt= / / aA<j> dxdt - V (f>(si, t)[a x ]i dt - / ^(s,, eft. (106) 
Jo Jo Jo Jo Jo i=1 ' Jo i=1 



Thus, using (|103[) . and noting that a is continuous at by construction, we arrive at 



(b(a)d t <t> + aA<p) dxdt + / b(a)(j)\ t= o dx = 
o Jo Jo 



(107) 



for all smooth (j>(x, t) such that (f>(x, T) = 0, as the weak formulation of SP n . 

This is identical to the problem considered in Chapter V.9 of [I], and existence and uniqueness 
of solutions follows by exactly the same argument used there, since one merely requires that b(a) be 
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piecewise smooth and monotonically increasing. We record this result as 

Lemma C.l For a given bounded, continuous initial datum ip(x), equation has a unique 

bounded solution o~(x,t). 
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